RECTANGULAR MIXED ELEMENTS FOR ELASTICITY WITH 
WEAKLY IMPOSED SYMMETRY CONDITION 



GERARD AWANOU 

Abstract. We present new rectangular mixed finite elements for linear elasticity. 
The approach is based on a modification of the Hcllingcr-Reissner functional in 
which the symmetry of the stress field is enforced weakly through the introduction 
of a Lagrange multiplier. The elements are analogues of the lowest order elements 
described in Arnold, Falk and Winther [ Mixed finite element methods for linear 
elasticity with weakly imposed symmetry. Mathematics of Computation 76 (2007), 
pp. 1699-1723]. Piecewise constants are used to approximate the displacement and 
the rotation. The first order BDM elements are used to approximate each row of 
the stress field. 



I. Introduction 

The theory of elasticity is used to predict the response of a material to applied forces. 
The unknowns in the equations are the stress field, a symmetric matrix field which 
encodes the internal forces and the displacement, a vector field. For various reasons, 
mixed finite elements where one approximates both the stress and displacement are 
the methods of choice. One seeks the stress in the space of symmetric matrix fields 
with components square integrable and with divergence, taken row-wise, also square 
integrable. The displacement is sought in the space of square integrable vector fields. 
The pair forms a unique saddle point of the Hellinger-Reissner functional. It is very 
difficult to construct at the discrete level, finite element spaces which satisfy Brezzi's 
stability conditions. These conditions provide sufficient conditions for the stability of 
mixed finite element methods. Indeed for several decades before the work of Arnold 
and Winter [10, 11] the existence of such elements was an open problem. These 
elements have been extended to rectangular meshes in two dimension [3, 17], three 
dimension [13] and on tetrahedral meshes [5, 1]. Despite their relative complexity, 
mixed finite elements with symmetric stress fields are useful in certain situations [25] . 
If one desires simpler elements, one is forced to turn to nonconforming elements. Non- 
conformity can be introduced by weakening the symmetry condition or by weakening 
the requirement that the stress field is L 2 integrable. We refer to [12] for a review on 
nonconforming elements with symmetric stress fields and other approaches to linear 
elasticity. 

Stable mixed finite elements with weakly imposed symmetry have been introduced 
in [2, 6, 26, 28, 27, 24, 7, 9, 15, 23, 22, 19], The purpose of this paper is to present 
elements with weakly imposed symmetry for rectangular meshes. Precisely, we will 
use piecewise constants to approximate the displacement and the rotation and 18 or 
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12 dimensional spaces to approximate the stress field. The simplest older element on 
rectangular meshes in two dimensions is the one of [24] with 11 degrees of freedom for 
the stress, piecewise constants to approximate the displacement but a 4 dimensional 
space to approximate the rotation. The advantage of our element is that the rotation 
can be eliminated by static condensation. In three dimensions as well, our elements 
are simpler than Morley's elements. 

The paper is organized as follows: after some preliminaries in the next section, we 
present our low order elements in two dimension and then in three dimension. We 
conclude with some remarks on higher order elements. 



Let Q be a simply connected polygonal domain of M n , n = 2,3, occupied by a linearly 
elastic body which is clamped on dQ. We denote as usual by L 2 (f2,R n ) the space of 
square integrable vector fields with values in R n and H k (K, X) the space of functions 
with domain K C R n , taking values in the finite dimensional space X, and with 
all derivatives of order at most k square integrable. We let if(div, f2, X) be the 
space of square-integrable fields taking values in X and which have square integrable 
divergence. For our purposes, X will be either M the space of n x n matrices, § 
the space of n x n symmetric matrices, R n , or R, and in the latter case, we simply 
write H k (X). The divergence operator is the usual divergence for vector fields which 
produces a matrix field when acting on a matrix field by taking the divergence of each 
row. We will also need the space if (curl, fl, R n ) of square-integrable fields with square 
integrable curl. We recall that in two dimension for a scalar function q, curl(g) = 
(d 2 q, —dig) and in three dimension 



For a vector field in two dimension or a matrix field in three dimension, the curl 
operator produces a matrix field by taking the curl of each row. The norms in 
H k (K, X) and H k (K) are denoted respectively by 1 1 • | \ H k and \ \- \\k- We use the usual 
notations of Vk{K,X) for the space of polynomials on K with values in X of total 
degree less than k and P^^-fT, X) for the space of polynomials of degree at most 
k\ in x and of degree at most k 2 in y. Similarly, P^^^^, X) denotes the space of 
polynomials of degree at most k\ in x, of degree at most k 2 in y and of degree at most 
k 3 in z. We write Vk, 'PkiM an d V^^m respectively when X — R. 
The solution (a, u) e if(div, f2, §) x L 2 (f2,M n ) of the elasticity problem can be char- 
acterized as the unique critical point of the Hellinger-Reissner functional 



The compliance tensor A = A(x) :§—?>§ is given, bounded and symmetric positive 
definite uniformly with respect to x G fl, and the body force / is also given. In the 
homogeneous and isotropic case, 



where I is the identity matrix and A and \i are the positive Lame constants. 



2. Preliminaries 



curl(<?i, q 2 , q 3 ) = (<9 2 g 3 - d 3 q 2 , -diq 3 + d 3 qi, diq 2 - d 2 qi). 
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To treat both two and three dimensional problems in a unified framework, one pos- 
sibility is to use finite element differential forms [8]. However, for n = 2, 3 a simple 
device will suffice. We define P to be M when n = 2 and P = IR 3 for n = 3. Then we 
define asr = Tyi — r 2 i for a 2 x 2 matrix and asr = (r 32 — r 23 , Ti 3 — r 3 i, r 2 i — Tyi)' 
in three dimension. For a symmetric matrix field, asr = 0. Next, we define EI to be 
M 2 when n = 2 and EI = M for n = 3. For the formulation with weakly imposed 
symmetry condition, a critical point of the extended functional 



J(a,v)+ [ 7] 
Jn 



asr 



is sought over ff(div,Q,M) x L 2 (Q,R n ) x f 2 (fi,P). The unique solution (a,u,i) 
satisfies 

(Aa, r) + (div r, u) + (as r, 7) = 0, r G if (div, tt, M) , 
(2.1) (div<7,u) = (/,«), t)6i 2 (fl,r), 

(as(T,g) = 0, ga 2 (fi,P). 

For the associated discrete system with finite element spaces Eft x 14 x Qh c if (div, f2, M) x 
L 2 (f2,IR n ) x L 2 (f2,P), the symmetry condition will be enforced only weakly. The 
Brezzi's conditions for stability are 

• There exists a positive constant c\ independent of h such that ||r||H(div) < 
Ci(At,t), if r G Eft, (divr, t>) = for all i> G 14 and (asr, g) = 0,Vg G <5h, 

• There exists a positive constant c 2 independent of /i such that V (v, q) G 14 x 
<5fc, 7^ (0,0), 3 r G Eft, r 7^ with (divr,w) + (asr,g) > c 2 ||r||# (div) (|H| L 2 + 
ll?IU 2 )- 

To fulfill these conditions, we construct Eft, 14 and such that 

1- div Eft C 14 

2- Given (v, q) G 14 x Q h , (t> , g) 7^ (0, 0), 3 r G Eft, r^0 such that 
(2-2) IMItf(div) < C{\\v\\ L 2 + \\q\\v), 

and divr = v, PQ h asr = q, where Pq h is the L 2 projection operator. 

The first Brezzi condition follows from the condition div Eft C 14 . It is easy to see 
that the second follows from condition (2) above. To construct elements which satisfy 
(1) and (2), we follow the constructive approach of Arnold, Falk and Winther, [7, 9], 
using discrete versions of the de Rham sequence. In addition to the spaces Eft, 14 and 
Q h , we also construct finite element spaces R h C if (div, f2, EI) and <d h C if (curl, f2, EI) 
in such a way that the following diagrams commute: 

if (div, ft, H) ^ L 2 {Q,F) — + 

Rh — > Qh — > 0, 
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if (curl, ft, H) ff (div, ft, M) L 2 {Q,R n ) 



r ^ curl „ div t r „ 

©ft — >■ — * 14 — ^ 0. 
We note that the commutativity of the far left side of the diagram above will not be 
used. For a finite dimensional space Xh, Hx h is a bounded projection operator. We 
recall that 

(2.3) U Xh v = v, V v G X h . 

Next, we define an operator S : C°°(ft,H) — > C°°(ft,H) which connects the two 
diagrams above. In two dimension, 5 is simply the identity operator, while in three 
dimension, for q = (%)ij=i,...,3, we define 

/ ?22 + ^33 -921 -931 

(2.4) S(q) = -qi2 911 + 933 -932 

\ -9i3 -923 911 + 922, 

In that case, S is also invertible with S(q) = tr(g)J — q T , 5'~ 1 (g) = l/2tr(g)f — q T , 
[15], where q T denotes the transpose of q, I is the 3x3 identity matrix and tr(g) 
denotes the trace of q. The following fundamental relation holds in both dimension: 

(2.5) ascurl q = — div S(q). 

We summarize the elements of the constructive approach of [7, 9] in the following 
theorem, the proof of which is reproduced below for convenience. 

Theorem 2.1. Under the commutativity assumptions 

(2.6) U Qh dwq = div U Rh q, Vg G C°°(ft,H), 

(2.7) divider = U Vh diva, Vcr G C°°(ft,M), 
and 

(2.8) n^^ne^-^n^, 

(2.9) ||n Sh M|| i2 < c||w|| H i, Vr G H\Q,M), 

(2.10) \\cmlU eh p\\ L 2 <c\\p\\ H i, VpGff^ft,!!). 
£/ie second Brezzi condition holds. 

Proof. By elliptic regularity, given t> G V^, 3 77 G ff^ft, M) such that 

(2.11) divi] = v and H^Hh 1 < IMU 2 - 
Given g G Q h C f 2 (ft,P), there exists G /^(f^H) such that 

(2.12) div0 = g - n Qh asH^r? and ||0||^i < C||g - Tl Qh as IIs^l | L 2. 
We set r = II 2h ?7 + curHle h S'~ 1 and by (2.7) and (2.3) we have 

divr = div Yl^ h r] = Hv h div 77 = Hv h v = v. 
By (2.5) and (2.6) it follows that 

Tl Qh as curl q = Tl Qh div Sq = div U Rh Sq, 
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We therefore have using (2.8), (2.6) and (2.3), 

U Qh asr = Yl Qh asIT Sh ?7 + Yl Qh as curl lie,, -S^V 
= U Qh asU^rj + div U^SUe^cf) 
= U Qh as n Sh ?? + div Tl Rh (f) 
= U Qh asU Eh r] + U Qh div0 
= U Qh as U^r] + U Qh q - U Qh as U Sh r] 
= q- 

It remains to prove the inequality (2.2). We have by (2.11) and (2.9) 

l|n Eh »7lU» <c|H|hi <c\\v\\ L 2, 

and by (2.11), (2.3), (2.11), (2.9) and (2.12) 

||curin ^-VlU 2 < \\S- l (f>\\m <C\\<j>\\ H i < \\q-U Qha sU^ h rj\\ L 2 

< C(\\q\\ L 2 + ||asn s ^|| L2 ) < C(\\q\\ L 2 + \\ V \\ H i) 

< C(\\q\\ L 2 + \\v\\ L 2). 

It follows that ||t|| L 2 = ||II Sfi ?7 + curin 0h 0|| L 2 < C(||g|| L 2 + ||t>|| L 2). Since divr = v, 
this proves the result. □ 

Let l~h denote a conforming partition of Q into rectangles of diameter bounded by h, 
which is quasi-uniform in the sense that the aspect ratio of the rectangles is bounded 
by a fixed constant. Let R = [0, l] n be the reference rectangle and let F : R — > R 
be an affine mapping onto R, F(x) = Bx + b, with b G R n and B a n x n diagonal 
matrix. Our goal in the next section is to construct spaces V h and 6^ such that 
the conditions of Theorem (2.1) hold. If (a, u,p) denotes the solution of problem (2.1) 
and (<Th,Uh,Ph) £ S k x 14 x 0i is the solution of the associated discrete system, the 
optimality condition 

, lk-<7fc||fT(div) + \ \u-U h \\ L 2 + 117-7/ilU 2 < C 'u&r h €Z h ,v h eV h ,p h eQ h 

(Ik - T h \\ H{div) + \\U - V h \\ L 2 + ||7 - p h \\ L 2), 

holds. 

As with [7, 5, 15], the following refined error estimates hold 

\W - cr h \\ H (div) + \\u h - Tl Vh u\\ L 2 + ||7 - 7 h || £a < C(\\a - U.^ h a\\ H (^y) + ||7 ~ n Q h 7lU 2 )> 
\\u - u h \\ L 2 < C(\\a - n Sh o-||^ (div) + ||7 - n Qh 7|| L 2 + ||u - U Vh u\\ L 2), 
|| div(er — (7^)11^2 = || div a" — Uv h divcr|| £ 2. 

3. TWO DIMENSIONAL ELEMENTS 

We recall the lowest order BDM element, 

(3.1) BDMi(K) = {q\q = p±(x,y) +rcurl(x 2 y) + s cnv\(xy 2 ),pi eV 1 xV 1 }, 

and an element q G BDMi(K) is uniquely determined by the conditions J e q ■ 
npi ds, for each edge e of K, Vpi G Vi(e). 
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We choose 14 = Vo(Th), Qh = T- > o(7h), with degrees of freedom the value at an interior 
point in each element K and 

Z K = {T,T(x,y) GM, (r^Taje^M^^U}. 

A matrix field r G is uniquely determined by the first two moments of rn on each 
edge, (2 x 2 x 4 = 16 degrees of freedom). The stress field space T, h is therefore the 
space of matrix fields which belong piecewise to T, K and have normal components 
which are continuous across mesh edges. 

We will also need the serendipity finite element space Sh, defined on a single element 
Kbj 

Sk = Vi{K) + span{aj 2 y,x|/ 2 }, 
and with degrees of freedom for q G Sk 

(1) the values of q at the vertices (4 degrees of freedom), 

(2) the average of q on each edge (4 degrees of freedom). 

It is not difficult to check that the sequence 

— >■ R A S K BDMi(K) ^ V (K) — >• 0. 

is exact. One checks that each space is mapped in the one that follows. Then one 
notes that the alternating sum of the dimensions is zero and that the polynomial de 
Rham sequence is exact. 

We therefore define the space as follows: on each element K, Qk = Sk x Sk 
and the space 6^ is the space of vector fields which belong piecewise to 6^ and are 
continuous across mesh edges. 

Finally we take for R h the lowest order Raviart-Thomas element, i.e. Rh = RT (Th)- 
We recall that RT (K) = V\ : o(K) x Vo,i(K) with degrees of freedom the average of 
the normal component of q G RT (K) on each edge. 

The projection operator n S/i is taken as the canonical interpolation operator and 
defined by 

J n Sh ((r)n ■ qds = J an ■ qds, for all edges e and for all q G V\(e) x Vi(e). 

Similarly we define H# h by 

J Uji h (q) ■ nds = J q ■ nds, for all edges e. 

It remains to define the interpolation operator Ile h . For this we first define 11^ : 
H\K,R 2 ) ->• Q K by 

H K ip(v) = for each vertex v of K, 

and Il£ : H\n,R 2 ) ->■ Q h by (n.° h T)\ K = U° k t. Next, let L h be a Clement interpola- 
tion operator [14, 18] which maps L 2 (f2,M) into 

{9 h e C°(ti)\e hlK eV^VK eT h }, 



and denote as well by the corresponding operator which maps L 2 (Q,WL 2 ) into the 
subspace Qh of continuous vector fields whose components are piecewise in Vi t i- We 
have 

(3.2) \\L h r - rl < ch m -mr\\ m , < j < 1, j < m < 2, 
with c independent of /i. We define our interpolation operator YIq h by 

(3.3) Ue h — n°(7 — L h ) + L h . 

Theorem 3.1. For the triple (X^, V^, O/J the conditions of Theorem (2.1) hold and 
we have the optimality condition (2.13). Moreover if a and u are sufficiently smooth, 

(3.4) \\a - a h \\H(div) + \\u - u h \\ L 2 + ||7 - 7 h || L 2 < C h\\u\\ 3 . 

Proof. Let q G C°°(f2, IR 2 ). We have using the definition of Uji h and Green's theorem, 
/ divilf^gcb = V] / div rfx = V] / U Rh q-nds 

Jn K JK K JdK 

= / q ■ nds = / div g, 

which proves (2.6). 

Next, let a - G C°°(f2,M). Again using the definition of and Green's theorem, 

/ diver — div H^ h a dx = / div(a — n Sh <7) <ir = / (cr — n 2;i 0")n rfs = 0, 

Jn K JK K JdK 

which proves (2.7). 

For q G C°°(f2,IR 2 ), put u = U^q. We have using the definition of 11° 

It follows that n flh (tx - g) = i.e. n flh II°9 = II flh g. Finally, n^n , = n^n°(7 - 
Lh) + n Rh L h = U Rh (I - L h ) + U Rh L h = U Rh , that is (2.8) holds. 

By the trace theorem, one shows that (Il Sh )|^ is bounded on i7 1 (/<", M). Moreover 
if we define for a matrix field M, P F {M)(x) = l/det(B)M(x)B T ,x = F(x), then it 
is not difficult to verify that P F ((U.^ h )\^a) = (Us h )\ K P F a, hence (2.9) follows from 
a standard scaling argument. 

Let p G H 1 ^, IR 2 ). We define its Piola transform by P F p = (PfPi, PfPi) where for 
a scalar function u, P F u — u o F -1 . 

Since curing/) G S^, 

| IcurHI^pl | L 2(f) <C / curing./) • hs 1 ds 

ecdit *=o Js 
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where e is an edge of dK. Next, curlg • n = dq/ds and using the definition of n^, 
Jcm:m\p-nds = J J^lfypds = 

J curing p ■ n sds = J-^-(Il ^p)sds = — J H°^pds = — J pds. 

By the trace theorem, it follows that 

||curin^p|| L2(f) < C\\p\\ h f, 
and scaling to an arbitrary rectangle K, we get 

|| cur\Il K p\\ L 2( K ) < C(/i _1 |p| , J ft: + C\p\ ljK )- 

We therefore have 

|| curHIe h p||L2 < || curin°(7 - L h )p\\ L 2 + 1 1 curl L h p\ | L 2 

< c^h-'l | (/ - L h )p\\ L , + 1 1 (/ - L h )p\\ H .) + c\\L hP \ \ H i 

< c IIpII//i, 

that is (2.10) holds. Since divS/j C Vh, the Brezzi conditions hold and the error 
estimates follow from the optimality error estimate from the theory of mixed methods, 
properties of the canonical interpolation operator for BDM elements, [16] p. 132, and 
error estimates of the L 2 projection operator. □ 

3.1. Simplified element of low order. Analogous to the simplified element of [7], 
we can develop elements simpler than the lowest order BDM type elements. The key 
point is that for (2.8) to hold, we only need 0^ to have normal components continuous 
across edges. We start the construction by taking as Q h the rectangular version of a 
space introduced by Fortin, [20] and [21] p. 153. The spaces Rh, Vh and Qh are the 
same. To define the space Qh, let i, j be the unit vectors in the x and y directions 
respectively. We put 

Pi = —x(l — x)(l —y)i 
P2 = -y(l - y)(l - x) j 
£>3 = x(l — x)yi 
p 4 = xy(l-y)j, 

and define on each element K, 

Qk = V lt i(K) x V hl {K) espan{p 1 ,p 2 ,p 3 ,p 4 } 
with degrees of freedom 

(1) the values of q at the vertices (4x2 = 8 degrees of freedom), 

(2) the average of q ■ n on each edge (4 degrees of freedom). 

The stress space Y, K is defined as 

(S'o(iO 7>oi(*o) ® SPan ^ CUrlpi ' curlp2 ' curl P3, curlp 4 }, 



where ( ^ (^1 ^> 0,1 l^-\ ) is the space of matrix fields with components in the i 
V'iM-frJ '0,1 \&)J 

„ . /x(l-x) (l-2x)(l-y)\ 

dicated spaces. Explicitly, we have curl pi = ( q q j,curlp 2 

\ /x(l-x) -(l-2x)y\ , 

', curlp 3 = I n n and curlp 4 



(-l + 2y)(l-x) - y (l- y )J>™ Fd ~ V 




x(l-2y) -y{l-y))' 

For r G ( ^ 1,0 /^-\ -t) ' 1 /^ ) > rn e ^o(e) x "P ( e ) on each edge e but (curlpj)n • t G 
Vi(e),i = 1, . . . , 4. The following degrees of freedom are unisolvent: 

(1) f rn ■ nds for each edge e 

(2) f e rn ■ tpds for each edge e and p G Vi(e). 

To see this, let r = rj + ai curlpi + a 2 curlp 2 + 03 curlp 3 + a 4 curlp4 G such that all 
the above degrees of freedom vanish. Since the normal component of (th, r i2 ),i = 1,2 
vanish on each edge, we have 



Since 



Tn = x(l - x)c iU T i2 = y(l - y)c i2 ,i = 1,2, c id G R,i,j = 1,2. 

T11 = V11 + aix(l -x) + a 3 x(l - x),rj u G V W (K) 
T12 = V12 + oi(l - 2x)(l -y)- a 3 (l - 2x)y,r} 12 G V i(K) 
T21 = V21 + fl2(-l + -x)- a 4 x(l - 2y),r} 21 G V W (K) 
T22 = V21 - a 4 y(l -y)- a 4 y(l - y),rj 22 G V i(K), 

we conclude that a\ = a 2 = a 3 = a 4 = and 77 = 0, that is: r = and the claim 
follows. 

From the approximation properties of the lowest order Raviart-Thomas element, the 
estimate (3.4) still holds. 

4. Three dimensional elements 

The de Rham complex in three dimensions is 

R ^ C°°(Q,R) ^4 C°°(Q,R 3 ) ^> C°°(Q,R 3 ) C°°{Q,R) — ■» 0. 
We choose the following form of BDM elememt, [16], p. 124 

GD.U.l/v): 'P;i A.R 3 ) + curlspan<| ( 1.1 I I I I I , I xz 2 I , I x 2 z 




y xy l j \x 2 yj \ J \ / \ / \ 
Clearly div BDM X {K) = V (K). We define = P (^) 3 and 

Y, K = { r, r(x, y, 2) G M, (r a , r i2 , r l3 ) G BDMi(K),i = 1, 2, 3 }. 
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The degrees of freedom on Vk are the values of each component at an interior point 
while a matrix field r in is uniquely determined by the moments of order and 
1 of rn on each face (3x3x6 degrees of freedom). 

We now define two spaces Sk and Uk such that the sequence below is exact. 

R A S K ^ U K ^\ BDMi(K) ^ V (K,R) — ► 0- 

The space Sk is not directly used in the construction but helped discover Uk- We 
take the space Sk as the three dimensional serendipity space of order 2 defined as 

Sk = Vi{K, K) + span{ x 2 y, x 2 z, xy 2 , xz 2 , y 2 z, yz 2 , xyz, x 2 yz, xy 2 z, xyz 2 }, 

with degrees of freedom 

(1) the values of q G Sk at the vertices (8 degrees of freedom), 

(2) the average of q G Sk on each edge (12 degrees of freedom). 

The unisolvency of these degrees of freedom is proven for example in [4]. We define 
the space Uk as 

Uk = V\^\{K, M 3 ) + span{ y 2 z, yz 2 , y 2 , z 2 } x span{ x 2 z, xz 2 , x 2 , z 2 } x span{ x 2 y, xy 2 , x 2 , y 2 }, 
with degrees of freedom for u G Uk, 

(1) the first two moments of u ■ t on each edge, where t is a tangential vector to 
the edge (12 x 2 = 24 degrees of freedom), 

(2) the average of u A n on each face with unit outward normal n (6 x 2 = 12 
degrees of freedom). 

It is not very difficult to verify that the sequence above is exact. One checks that 
each space is mapped in the one that follows. Then one notes that the alternating 
sum of the dimensions is zero and that the polynomial de Rham sequence is exact. 
We then only need to verify either that the kernel of the curl operator is the image 
of the grad operator or that the kernel of the div operator is the image of the curl 
operator. We verify the last one. Let u G BDMi(K) such that div-u = 0. We 
write u = w + curlew G V\[K, M 3 ) and z in the span of the extra monomials in 
the definition of BDMi(K). Note that z G Uk and div-u = divw = 0. By the 
exactness of the polynomial de Rham sequence, w = curia, a G V2{K,M. 3 ). Since for 
a, (3, 7 G K, curl(o;x 2 , (3y 2 , ^z 2 ) = 0, we may assume that a G Uk which completes the 
proof of the claim. 

We can now describe the space Q h as 

®h = {q, q(x, y, z) G M, (q a , q i2 , q i3 ) EU h ,i = 1, 2, 3 }, 
with the degrees of freedom for q E Oh 

(1) f qts l ,i = 0, 1 for each edge e, where t is a tangential vector to the edge 
(12 x 2 x 3 = 72 degrees of freedom), 

(2) Jj.qAndxf for each face / with unit outward normal n (6 x 2 x 3 = 36 degrees 
of freedom). For a matrix field q with row vectors a;, i — 1, 2, 3, qAn is defined 
as the matrix field with rows qi A n, i — 1, 2, 3. 
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Next we define the space Qh- We take Qk = 'Po(K) 3 with degrees of freedom the 
values of each component at an interior point. 

Finally we describe the space Rh as 

{ q, q(x, y, z) G M, (q iU q i2 , q a )\ K G RT (K),i = 1, 2, 3 }, 

where 

RT (K) = V 1>0 , {K) x V , 1>0 (K) x P , ,i(K), 

is the lowest order Raviart-Thomas element in three dimensions with degrees of free- 
dom the average of the normal component on each face, (1 x 1 x 6=6 degrees of 
freedom) . 

4.0.1. Unisolvency. The unisolvency of the degrees of freedom for Vk, and Sk 
are well known. Similarly unisolvency for the degrees of freedom of Rh is immediate. 
We only study the case of Uk- Let v = (v±, v 2 , v 3 ) G Uk and assume that all degrees 
of freedom vanish. We show that v\ = 0. On each edge e, v - t G V\{e) and hence we 
get v ■ t — on each edge. This implies that on the face z = for example, 

V! = y(l - y)w 1 ,w 1 G Via 
v 2 = x(l - x)w 2 , w 2 G Vo,i- 

However, if w\ has a linear term in x, xy 2 would be the highest degree monomial in 
vi. We conclude that wi is constant. The face degrees of freedom imply that the 
average of w\ vanish on the face z — 0, that is: w 1 = 0. Similarly w 2 = 0. We 
conclude that v has expression 

vi = 2/(1 - y)z{l - z)r u 
v 2 = x(l - x)z(l - z)r 2 , 
v 3 = x(l - x)y(l - y)r 3 , 

for constants r 1 ,r 2 and r 3 which must vanish given the form of the highest degree 
monomial in the expression of v^i = 1, 2, 3. 

4.0.2. Definition of interpolation operators. For q G C°°(f2, M), we define Il^ h by 

J (UR h q)ndx = J qndx, for all faces /. 
The interpolation operator Us h is defined by 

/ n Sh (o-)n • qds = an - qds, for all faces / and for all q G Vi(f) x V\(f) x Vi(f). 
Jf Jf 

It remains to define the interpolation operator He h - For this we first define 11^ : 
H l (K, M) — >■ Q K by 

/ (U° K q) t s i ds = 0, % = 0, 1 for each edge e C dK, 
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and Tl° h : H\tt, M) ->■ O h by (Tl° h T)\ K = TI° k t. Next, let L h be a Clement interpolation 
operator [14, 18] which maps L 2 (f2,M) into 

{ e h g c°(Q) | e h \ K g Vi,i,i,vk g r h }, 

and denote as well by the corresponding operator which maps L 2 (f2,M) into the 
subspace of 0^ of continuous matrix fields whose components are piecewise in Pi,i,i- 
We have 

(4.1) \\L h r - rl < ch m -mr\\ m , < j < 1, J < m < 2, 
with c independent of h. We define our interpolation operator He^ by 

(4.2) U.Q h — n°(7 — L h ) + L h . 

4.0.3. Commutativity and surjectivity assumptions. The commutativity assumption 
(2.6) and (2.7) are proven as in the 2D case. We verify the surjectivity assumption 
U Rh SUe h = Tl Rh S. We first show that IL Rh SUe h = U Rh S. For this let q G C°°(fi, M), 
put u = q — 11° 9. We need to show that U Rh Su> = 0, that is 



Since U° h w = 0, 



J {Su){x)nalxf = 0, for each face /. 



u A n = 0, for each face /. 



/ 



Next for q = (9^)^=1,2,3, 

(913^1 - 9ll% 9ll% - 912^1 9l2% - 9l3%\ 
923% - 921% 921% - 922% 9 22 % ~ 923% , 
933% - 93in 3 931^2 - 9 3 2% 9 32 % - 933% / 

and 

( 922% + 933% - 921% - 931% \ /~(q A n) 22 + (q A n) 31 
~9i2% + 911% + 933% - 932% = (g A n)i2 - (g A n) 3S 
-9i3% - 923% + 9n% + 922% / \-(g A n) u + (9 A 71)23, 

This shows that J f uAn = implies J^(Su)n = and the result follows using the 
definition of 11^. 

We notice that for 9 G 6/j, for the surjectivity assumption to hold, the following de- 
grees of freedom were not used: f 912%— 9i3% dxf = Jj(qAn) 13 , J 923%— 921% dxf — 
/ (9 A 71)12, Jj 931% — 932% dxf = / (9 A n) 32 . However since the faces of a rectangle 
are parallel to the axes, one of these degrees of freedom is identically zero for each 
face, hence two degrees of freedom per face are unnecessary. 

4.0.4. Boundedness of the interpolation operators. By the trace theorem, one shows 
that (IIs h )|^- is bounded on H l (K, M). Moreover if we define for a matrix field 
M, P F (M)(x) = l/det(B)M(x)B T ,x = F(x), then it is not difficult to verify that 
^((IlEjJI^-a) = (H-x^IkPfO'i hence (2.9) follows from a standard scaling argument. 

Let p G H l (K, M?). We define its Piola transform by Ppp = (PfPi, PfP2, PfPz) where 
for a scalar function u, Pfu = u o F^ 1 . 
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Since curlli^p gE^-, 

| Icurlll^pl | jL 2(f.) 

fCdK i=0 * 

where / is a face of dK. Next, using the definition of U ^, for q G V\ i(f) x V\ i(/) x 

(curl(n^.p)n) -qdxf = J(Il ^p) AhVqdxj = J p A hVq dx /. 

By the trace theorem, it follows that 

||curin^p|| L2(f) < C||p|| ljf , 
and scaling to an arbitrary rectangle K, we get 

|| curing | L 2 (A -) < C(/i _1 |p|o,^ + C\p\i )K )- 

We therefore have 

|| curin eh p||i2 < || curing/ - L h )p\\ L 2 + 1 1 curl L h p\ \ L 2 

< cihT^ | (/ - L h )p\\ L 2 + 1 1 (/ - L h )p\\ m ) + c\\L h p\ \ m 
<c\\p\\m, 



that is (2.10) holds. Since div£ fe C Vh, the Brezzi conditions hold. From the op- 
timally error estimate from the theory of mixed methods (2.13), properties of the 
canonical interpolation operator for BDM elements, [16] p. 132, and error estimates 
of the L? projection operator, we have the following error estimate. 

Theorem 4.1. For the triple (E^, 14, O^) the conditions of Theorem (2.1) hold and 
we have the optimality condition (2.13). Moreover if a and u are sufficiently smooth, 

(4.3) \\a - cr h \\ Hidiv) + \\u - u h \\ L 2 + ||7 - ^ h \\ L 2 < C h\\u\\ 3 . 

5. Higher order elements 

Except the simplified element in two dimension, the elements we have described do 
not have optimal rate of convergence for the stress. It does not seem possible to 
simplify the three dimensional element using the framework described here. In two 
dimension, for higher order approximation, if (div) elements can be constructed based 
on the sequence, 

— > R Vk+i,k+i — > Vk+i,k x Vk,k+i — > Pk,k — > 0. 
Take Vh to be the space of piecewise continuous vector fields which belong locally 
to Pk t k{K) x Pk t k{K), Qh the space of piecewise continuous functions which belong 
locally to Q K = P k -i,k-i(K) and S x = { r G M, (r a , r i2 ) G V k +i,k x V k ,k+u * = 1> 2 } 
with degrees of freedom 



(1) f rn ■ p k ds, for each edge e of K, V p k G Vk{e), 
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for k > 1. The space Rh is taken to be the Raviart-Thomas space of order k — 1 
and finally the space Oh is the space of continuous vector fields with components in 
Vk+i,k+i{K) on each element K. Again, there one does not have optimal convergence 
rate for the stress. We leave the details of the three dimensional analogue to the 
interested reader. 
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